program main
  use precision,only:p_
  use constants,only: m,pi   ! m is the number of velocity grid number.
  implicit none
  integer:: i,j,nstep,flag
  real(p_):: dv,dt,zi
  real(p_):: v(m),q0(m),q(m),fm(m)
  real(p_):: aa(m),bb(m),cc(m),col(m),dvv(m),djj(m)
  real(p_):: a(m),b(m),c(m),r(m)
  real(p_):: cond     !electrical conductivity
  real(p_):: tmp,qp(m),c10(m),sum
  real(p_):: q_old(m),det_q(m),det
  real(p_),parameter:: epsilon=1.0e-6
  ! this namelist provides the parameter Zi and the step number in time.
  namelist /input/ nstep,zi  
  open(40,file='input_file')
  read(40,input)
  !create grid and give maxwellian distribution and initial distribution for q.
  dv=15.0_p_/(m-1)
  do j=1,m
     v(j)=0.0_p_+(j-1)*dv
     fm(j)=1./sqrt((2.*pi)**3)*exp(-0.5*v(j)**2)
     !fm(j)=1./sqrt((pi)**3)*exp(-v(j)**2)
     q0(j)=1.0 
  enddo

  !calculate the difusion coefficient for maxwellian background.
  call get_d_f(dv,v,fm,dvv,djj) 
  ! calculate the coefficient of the differential equation (adjoint equation)
  call get_coefficient(zi,dv,v,fm,dvv,djj,aa,bb,cc) 

  !time step size. Larger size can be used thanks for the implicit scheme.
  dt=1000.0_p_   
  q=q0 ! initial condition for q
  open(30,file='spitzer.txt')  ! this file stores the steady distribution function.
  flag=1
  ! advance in time
  do i=1,nstep
     ! calculate the collision term between Maxwellian and first Legendre harmonics.
     call get_col(dv,v,fm,q,col)
     ! constuct the tridiagonal matrix
     call get_matrix(dt,dv,v,q,aa,bb,cc,col,a,b,c,r)
          q_old=q
     ! sovle the tridiagonal matrix equation
     call tridag(a,b,c,r,q,m)
      det_q=q-q_old
     det=0._p_
     do j=1,m
        det=det+abs(det_q(j))
     enddo
     if(det.le.epsilon)  then
        flag=0
        exit
    endif    
  enddo
 if(flag==1)   write(*,*) '****Not arrive at the desired precision*****'
 write(*,*) 'iteration number=',i, 'det=',det
  !output the steady distribution function
  ! gp is the analytic spitzer function.
  do j=1,m
     tmp=v(j)/sqrt(2.0_p_)
     qp(j)=-tmp**2*(0.6_p_+1.41_p_*tmp-0.66_p_*tmp**2+0.134_p_*tmp**3)
  enddo
  do j=1,m
     write(30,100) v(j),q(j),4.*qp(j),v(j)/sqrt(2.0)
  enddo
100 format(4f20.6)
  close(30)
  ! calculate the electrical conductivity
  call get_conductivity(dv,v,fm,q,cond)
  write(*,*) 'conductivity=', cond
  !the tabulated spitzer function differs a factor 4 from spitzer function defined here.
  !q=4._p_*qp  
  ! calculate the electrical conductivity
  !call get_conductivity(dv,v,fm,q,cond)
  !write(*,*) 'conductivity for z=1, and use analytic spitzer function', cond

  !prove that C(f1,fm)+C(fm,f1) conserves momentum.
  q=1.*qp
  call get_col(dv,v,fm,q,col)
  call get_c10(dv,v,q,fm,c10)
  sum=0.0  
  do j=2,m-1
     sum=sum+(col(j)+c10(j))*fm(j)*v(j)**3*dv
  enddo
  write(*,*) 'sum=',sum

end program main



subroutine get_d_f(dv,v,fb,dvv,djj)
  ! This routine calculates the difusion coefficient,dvv,djj, for isotropic background, fb.
  ! dvv:  difusion coefficient Dvv
  ! djj:  difusion coefficient D_theta_theta
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: dv,v(m),fb(m)
  real(p_),intent(out)::dvv(m),djj(m)
  integer::i,j
  real(p_):: sum1,sum2

  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+v(j)**4*fb(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+v(j)*fb(j)*dv
     enddo
     dvv(i)=4.*pi/3.*(sum1/v(i)**3+sum2)
  enddo

  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+v(j)**2/(2.*v(i)**3)*(3.*v(i)**2-v(j)**2)*fb(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+v(j)*fb(j)*dv
     enddo
     djj(i)=4.*pi/3*(sum1+sum2)
  enddo
end subroutine get_d_f

subroutine get_coefficient(zi,dv,v,fm,dvv,djj,aa,bb,cc)
  !This routine calculate the coefficient of the adjoint equation.
  !aa(m) for coefficient before second derivative
  !bb(m) for coefficient before first derivative
  !cc(m) for coefficient before zero derivative
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: zi,dv,v(m),fm(m),dvv(m),djj(m)
  real(p_),intent(out)::aa(m),bb(m),cc(m)
  integer::i,j
  real(p_):: sum1,sum2,tmp

  do i=2,m-1
     aa(i)=dvv(i)
  enddo
  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+v(j)**4*fm(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+v(j)*fm(j)*dv
     enddo
     tmp=4*pi/3.*(-1./v(i)**4*sum1+2./v(i)*sum2) 
     bb(i)=tmp-v(i)*aa(i)
  enddo

  do i=2,m-1
     cc(i)=-2.*djj(i)/v(i)**2-zi/v(i)**3
  enddo

end subroutine get_coefficient

subroutine get_col(dv,v,fm,q,col)
  !This routine calculates the collision term between maxwellian and first Legendre harmonics.
  !This is an implementation of Eq.(34) in Karney's paper [Comp. Phys. Rep. 4(3-4),183-244, 1986].
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: v(m),fm(m),q(m),dv
  real(p_),intent(out)::col(m)
  integer::i,j
  real(p_):: sum1,sum2,tmp

  do i=2,m-1
     sum1=0.
     do j=2,i  !Use rectangle integration formula for the velocity integration
        tmp=v(j)**3/(5.*v(i)**2)-v(j)/(3.*v(i)**2)
        sum1=sum1+v(j)**2*tmp*fm(j)*q(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        tmp=v(i)**3/(5.*v(j)**2)-v(i)/(3.*v(j)**2)
        sum2=sum2+v(j)**2*tmp*fm(j)*q(j)*dv
     enddo
     col(i)=4.*pi*(fm(i)*q(i)+sum1+sum2)
  enddo
end subroutine get_col

subroutine get_matrix(dt,dv,v,q,aa,bb,cc,col,a,b,c,r)
  !This routine is to construct the tridiagonal matrix.
  use precision,only:p_
  use constants,only: m
  implicit none
  real(p_),intent(in):: dt,dv,v(m),q(m),aa(m),bb(m),cc(m),col(m)
  real(p_),intent(out):: a(m),b(m),c(m),r(m)
  integer:: i

  ! the following loop is to calculate 3 diagonal lines in the tridiagonal matrix
  do i=2,m-1
     a(i)=-aa(i)*dt/dv**2+bb(i)*dt/(2.*dv)
     b(i)=1.+2*aa(i)*dt/dv**2-cc(i)*dt
     c(i)=-aa(i)*dt/dv**2-bb(i)*dt/(2.*dv)
     r(i)=q(i)+dt*col(i)-dt*v(i)    

  enddo
  ! boundary conditions:
  ! the following 3 lines are to implement q(1)=0
  b(1)=1._p_
  c(1)=0._p_
  r(1)=0._p_
  ! the following 3 lines are to implement q''(m-1)=0
  a(m)=-b(m-1)-2*a(m-1)
  b(m)=a(m-1)-c(m-1)
  r(m)=-r(m-1)
end subroutine get_matrix


SUBROUTINE tridag(a,b,c,r,u,n)
  !This routine is to solve the tridiagonal matrix. This routine is from Numerical Recipes, 
  !with some minor modification.
  !(C) Copr. 1986-92 Numerical Recipes Software ,4-#.
  use precision,only:p_
  implicit none
  !  integer,PARAMETER:: NMAX=5000
  INTEGER,intent(in):: n
  REAL(p_),intent(in):: a(n),b(n),c(n),r(n)
  real(p_),intent(out)::u(n)
  INTEGER j
  ! REAL(p_) bet,gam(NMAX)
  REAL(p_) bet,gam(n)
  if(b(1).eq.0.) stop 'tridag: The tridiangle matrix equation need rewriting'
  bet=b(1)
  u(1)=r(1)/bet
  do  j=2,n
     gam(j)=c(j-1)/bet
     bet=b(j)-a(j)*gam(j)
     if(bet.eq.0.) stop 'tridag failed'
     u(j)=(r(j)-a(j)*u(j-1))/bet
  enddo
  do  j=n-1,1,-1
     u(j)=u(j)-gam(j+1)*u(j+1)
  enddo
END SUBROUTINE tridag


subroutine get_conductivity(dv,v,fm,q,cond)
  !This routine calculates the electrical conductivity from the perturbed distribution function.
  ! perturbed distribution function is q*fm
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: dv,v(m),fm(m),q(m)
  real(p_),intent(out):: cond
  integer:: j
  real(p_):: sum

  sum=0._p_
  do j=2,m
     sum=sum+v(j)**3*fm(j)*q(j)*dv
  enddo
  cond=-4*pi/3*sum
end subroutine get_conductivity


subroutine get_c10(dv,v,q,fm,c10)
  !This routine calculates the collision term of first Legendre harmonics off Maxwellian
  !backgound, c10=C(f1,f0)/[fm*cos(theta)]*t0, where t0=1/nu0, nu0 is the thermal collision frequency
  use precision,only:p_
  use constants,only: m,pi
  implicit none
  real(p_),intent(in):: dv,v(m),q(m),fm(m)
  real(p_),intent(out)::c10(m)
  integer::i,j
  real(p_):: sum1,sum2,tmp,dvv(m),djj(m),bb(m),qp(m),qpp(m)
  !calculate the difusion coefficient for maxwellian background.
  call get_d_f(dv,v,fm,dvv,djj) 

  do i=2,m-1
     sum1=0.
     do j=2,i
        sum1=sum1+v(j)**4*fm(j)*dv
     enddo
     sum2=0.
     do j=i+1,m
        sum2=sum2+v(j)*fm(j)*dv
     enddo
     tmp=4*pi/3.*(-1./v(i)**4*sum1+2./v(i)*sum2) 
     bb(i)=tmp-v(i)*dvv(i)
  enddo

  do j=2,m-1
     qp(j)=(q(j+1)-q(j))/dv ! first derivative
     qpp(j)=(q(j-1)-2*q(j)+q(j+1))/dv**2 ! first derivative
     c10(j)=dvv(j)*qpp(j)+bb(j)*qp(j)-2.*djj(j)/v(j)**2*q(j) !C10=C(f1,fm)
  enddo

end subroutine get_c10
